TD3 - 17/05/21 - Simulation et Monte Carlo

  • Exercice 7 : MCMC
  • Exercice 9 : Cross-Entropy Method

Exercice 7 : MCMC

Rappel :

image.png

7.1.

  • Quelle loi reconnaissez-vous ?

image.png

7.2.

  • Calculer les lois conditionnelles de chaque composante.

image.png

image.png

  • Mettre en oeuvre le Gibbs sampler correspondant

image.png

image.png

Code

In [3]:
from scipy import stats
import numpy as np

def rand_gauss_gibbs(n, rho):
    
    # initialisation : nous pouvons choisir autre valeur
    X = np.zeros((n, 2))
    
    # itération gibbs sampler
    for t in range(n-1):
        X[t, 0] = stats.norm.rvs(rho*X[t, 1], np.sqrt(1-rho**2), 1)
        X[t+1, 1] = stats.norm.rvs(rho*X[t, 0], np.sqrt(1-rho**2), 1)
    
    # ne pas oublier la dernière simulation de X1
    X[n-1, 0] = stats.norm.rvs(rho*X[n-1, 1], np.sqrt(1-rho**2), 1)
    
    return(X)

rn = rand_gauss_gibbs(10, 0.6)
In [4]:
rn
Out[4]:
array([[-0.20305043,  0.        ],
       [ 0.61230731,  1.28002217],
       [-0.92959718,  0.61371824],
       [-0.71798956, -1.72158907],
       [-0.21810525, -0.74227182],
       [-2.23561418, -0.99432847],
       [-0.26211769,  1.19139452],
       [ 0.78155815, -0.38278824],
       [ 1.79929744,  1.01073237],
       [ 0.21971534,  0.78874589]])
  • Représenter dans le plan : des contours de la densité cible et l'évolution de la chaîne de Markov simulée
In [5]:
import matplotlib.pyplot as plt

# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_gibbs_contour(n, rho):
    
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    
    # contours de la densité cible
    rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
    
    x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
    pos = np.dstack((x, y))
    
    plt.contour(x, y, rv.pdf(pos))
    plt.title(f"rho = {rho}")
    plt.xlabel(f"X1")
    plt.ylabel(f"X2")

    # représentation de la chaine de Markov simulée
    plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
    return None
In [6]:
plot_gibbs_contour(n = 500, rho = 0.5)
  • Que se passe-t-il quand $\rho \to 1$ ?
In [7]:
n = 100
rhos = [0, 0.2, 0.5, 0.7, 0.9, 0.99]

plt.figure(figsize = (16, 10))

for i, rho in enumerate(rhos):
    plt.subplot(2, 3, i+1)
    plot_gibbs_contour(n = n, rho = rho)

lorsque $\rho \to 1$, $X_2 = X_1$ p.s.

rq : on peut montrer (et voir) facilement que lorsque $\rho \to -1$, $X_2 = -X_1$ p.s.

Nous utilisons les package python statsmodels pour évaluer la performance MCMC

In [8]:
from statsmodels.graphics.tsaplots import plot_acf
import pandas as pd

n = 5000

def plot_gibbs_ACF(n, rho, dim = 0):
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    # plot
    plot_acf(rn[:,dim])
    plt.title(f"Autocorrelation pour X{dim+1}")
    return None

def plot_gibbs_trace(n, rho, dim = 0):
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    # plot
    plt.plot(np.arange(n), rn[:,dim])
    plt.title(f"Trace X{dim+1}")
    return None
In [9]:
rho = 0.5
plot_gibbs_ACF(n, rho, dim = 0)
plot_gibbs_ACF(n, rho, dim = 1)
plt.show()
In [10]:
rho = 0.5

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()

Que se passe-t-il si $\rho$ proche de 1 ?

In [11]:
plot_gibbs_ACF(n, 0.99, dim = 0)
plot_gibbs_ACF(n, 0.99, dim = 1)
plt.show()
In [12]:
rho = 0.99

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()

7.3.

  • Programmer un algorithme RWHM (Hastings-Metropolis avec marche aléatoire)
In [13]:
def RWMH(n, rho, sigma):
    X = np.zeros((n,2))

    for t in range(n-1):
        
        # proposition
        x = stats.multivariate_normal(X[t,:], sigma).rvs(1)
        
        rvn = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
        
        if np.log(np.random.rand()) < np.log(rvn.pdf(x)) - np.log(rvn.pdf(X[t,:])):
            X[t+1,:] = x
        
        else:
            X[t+1,:] = X[t,:]

    return X

rho = 0.8
tau = 0.8
X = RWMH(n, rho, tau * np.eye(2))
X[:10]
Out[13]:
array([[ 0.        ,  0.        ],
       [-1.40641682, -0.34943386],
       [-1.40641682, -0.34943386],
       [-1.40641682, -0.34943386],
       [-2.31520794, -1.35749994],
       [-2.29135821, -1.55224073],
       [-2.29135821, -1.55224073],
       [-2.29135821, -1.55224073],
       [-0.30236218, -0.87273014],
       [-0.30236218, -0.87273014]])
In [37]:
# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_RWMH_contour(n, rho, tau, sigma = np.eye(2)):
    
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    
    # contours de la densité cible
    rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
    
    x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
    pos = np.dstack((x, y))
    
    plt.contour(x, y, rv.pdf(pos))
    plt.title(f"rho = {rho}, tau = {tau}")
    plt.xlabel(f"X1")
    plt.ylabel(f"X2")

    # représentation de la chaine de Markov simulée
    plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
    return None

def plot_RWMH_ACF(n, rho, tau, sigma = np.eye(2), dim = 0):
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    # plot
    plot_acf(rn[:,dim])
    plt.title(f"Autocorrelation pour X{dim+1} avec rho = {rho}, tau = {tau}")
    return None

def plot_RWMH_trace(n, rho, tau, sigma = np.eye(2), dim = 0):
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    # plot
    plt.plot(np.arange(len(rn)), rn[:,dim])
    plt.title(f"Trace X{dim+1} avec rho = {rho}, tau = {tau}")
    return None
  • voici ci-dessous les paramètres que nous allons analyser :
In [15]:
taus = [0.001, 1, 5, 100]
rhos = [0.2, 0.9, 0.99]
In [16]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = n, rho = rhos[0], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
In [17]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = n, rho = rhos[1], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
In [18]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = n, rho = rhos[2], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
In [19]:
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
In [20]:
plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()

Ne tenons plus compte des doublons :

In [34]:
def RWMH(n, rho, sigma):
    X = [[0, 0]]

    for t in range(n-1):
        
        # proposition
        x = stats.multivariate_normal(np.array(X[-1]), sigma).rvs(1)
        
        rvn = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
        
        if np.log(np.random.rand()) < np.log(rvn.pdf(x)) - np.log(rvn.pdf(X[-1])):
            X.append(x)

    return np.array(X)
In [39]:
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log
  # This is added back by InteractiveShellApp.init_path()
In [41]:
plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()

Bilan des courses :

  • Gibbs sampling est un bon outil pour construire des chaînes de Markov

  • Cependant lorsque les v.a marginales sont très corrélées ($\rho \to 1$) nous avons des anomalies : loi simulée différente de la loi cible

  • nous pouvons améliorer la performance de la simulation en utilisant un algorithme RWHM (où $\Sigma$ de la loi auxilière $\mathcal{N}(\mu, \Sigma)$ doit être bien choisie, pas "trop grand", pas "trop petit").

Exercice 9 : Cross-Entropy

image.png

image.png

image.png

image.png

image.png

Code :

In [ ]:
from numpy.random import multivariate_normal

def rosenbrock(x):
    y = 0
    for i in range(len(x)-1):
        y = y + 100 * (x[i+1]-x[i])**2 + (x[i]-1)**2
    return(y)

def optCE(fun = rosenbrock, n = 1000, d = 3, eps = 1e-8, max_iter = 1000, tau = 0.01):
    """
    this function compute the arg max of the function 'fun' by cross-entropy method
    """
    mu = np.random.rand(d)
    sigma = 10*np.eye(d)
    
    t = 0
    while True:
        t += 1
        
        # step 1 : sample Yi
        Y = multivariate_normal(mean = mu, cov = sigma, size=n)

        # step 2 : pick the top Yi that minimize the function
        samples = - np.array(list(map(fun, Y)))
        n_top = round(tau * n)
        top_ind = np.argsort(samples)[::-1][:n_top]
        Y_top = Y[top_ind, :]

        # step 3 : estimate by MLE mu and sigma
        mu = np.mean(Y_top, 0)
        sigma = np.cov(Y_top.T)
            
        # step 4 : stopping criter
        if np.max(sigma) < eps or t >= max_iter:
            print(t)
            print(mu)
            break
    return mu
In [ ]:
optCE(n = 1000, d = 3)
In [ ]:
optCE(n = 1000, d = 5)
In [ ]:
optCE(n = 1000, d = 10)
In [ ]:
optCE(n = 10000, d = 10)

Bilan des courses :

  • la méthode CE peut-être utilisé soit pour une estimation soit pour une optimisation
  • pour une optimisation, on doit tenir en compte plusieurs paramètres :

    • si $n$ (taille de l'échantillon) devient grand alors l'optimisation est plus précise
    • si $d$ (dimension) devient grand alors l'optimisation est plus compliqué
  • l'initialisation de $\theta_0$ est bien sûr très important, ici $\theta_0 = (\mu_0, \sigma^2_0)$ et :

    • idéalement $\mu_0$ doit être proche de la vraie valeur (sinon on choisit une valeur déterministe ou aléatoire)
    • $\sigma_0$ doit être assez grand afin de considérer toutes les possibilités